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Neutron star crusts are studied using a classical molecular dynamics model developed for heavy 
ion reactions. After the model is shown to produce a plethora of the so-called "pasta" shapes, 
a series of techniques borrowed from nuclear physics, condensed matter physics and topology are 
used to craft a method that can be used to characterize the shape of the pasta structures in an 
unequivocal way. 



PACS numbers: PACS 24.10.Lx, 02.70.Ns 



I. INTRODUCTION 

Neutron stars created in the death of a massive star 
are composed of a dense core containing an excess of 
neutrons over protons, thus justifying the name. With 
a mass between 1 and 3 solar masses and a radius of 
about 10 km, the stars are topped with a crust of about a 
kilometer thick where the /3 decayed-produced neutrons 
form neutron-rich nuclear matter immersed in a sea of 
electrons. The crust density ranges from normal nuclear 
density (~ 3 x 10 14 g/cm 3 ) at a depth of about 1 km, 
to the neutron drip density (~4x 10 11 g/cm 3 ) at about 
1/2 km, to a light mix of neutron-rich nuclei with densi- 
ties decreasing practically to zero in the neutron star en- 
velope. Likewise the proton-to-nucleon ratio also varies 
from ~ 0.25 to ~ 0.5 through the crust [l[, and the tem- 
perature from cold nuclear matter to about 1 MeV. The 
study of the structure of such crust, is the purpose of the 
present work. 

Studies of low density nuclear matter have found that 
the attractive- repulsive interplay of nuclear and Coulomb 
forces drive low-density nuclear matter to form non- 
uniform structures which are collectively known as "nu- 
clear pasta" . Such arrangements go from condensed 
phases with voids filled with nuclear gas, to "lasagna- 
like" layers of nuclei surrounded by gas, to "spaghetti- 
like" rods embedded in a nuclear gas, to ever decreasing 
"meatball-like" clumps which practically dissolve into a 
gaseous phase 0. 

Early investigations have used static models which rely 
mostly on energy considerations to determine the struc- 
tures that are most energetically favorable. Among the 
various theories proposed, the ones used most recurrently 
arc the Compressible Liquid Drop Model @-@], the Ex- 
tended Thomas- Fermi Model 0, [a-Q , and the Hartree- 
Fock method 

On the other hand, there are dynamical models that 
go beyond mean fields to predict the formation of the 
pasta phases as an asymptotic equilibrium state resulting 
from an evolution of a dynamical system. The most used 
methods are the semiclassical molecular dynamics (T3 - 



[lU and the quantum molecular dynamics p7l - l22l |. 

On the general composition of the pasta, most models 
agree on the formation of varying structures at subnormal 
densities but not on how the sequence the phases arises. 
Since the physical mechanism responsible for the phase 
transition pattern is a subtle interplay between Coulomb 
and nuclear energies which varies only a few keV/ fm 3 
between phases, the precise transition pattern is easily 
altered by the ingredients of the theoretical models. The 
fact that the Coulomb interaction between the electron 
sea and the nucleons -whose screening effect stabilizes the 
overall system modifying its structure- must be treated 
under different approximations in different models [22[ 
complicates any cross-model comparison even more. 

An additional problem which makes comparisons be- 
tween models difficult or impossible in some cases, is 
the lack of a quantifiable characterization of different 
pasta phases. The identification of phases has been done 
mainly through visual inspections of snapshots of spa- 
tial nucleon distributions obtained from calculations [23| . 
The QMD approach of reference [llj, for example, pro- 
duces nuclear holes, slabs, cylinders and spheres similar 
to those predicted by the Thomas Fermi model @, @, 01 , 
but at different densities and temperatures. 

In spite of this, the pasta phases have been character- 
ized globally. For instance, static models have been used 
to calculate average densities [|| and volume fractions 
of the different phases 1241 . Pasta bulk properties, such 
as the shear viscosity [25j and diffusion coefficients (26j 
have also been obtained using molecular dynamics sim- 
ulations. Refined studies have used radial correlation 
functions to characterize the nucleon distributions [lj| 
and the pasta structure factor to study charge density 
fluctuations [27^ 

More recently, shape characterizations were attempted 
both in dynamical simulations and with static mod- 
els; the former use topological measures such as 
the Minkowski functionals and the Euler characteris- 
tics [23|, while the latter modified the liquid drop model 
with a curvature correction to detect structure shape 
changes p8j . 

Thus the motivation of the present study: how to 
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achieve a precise enough characterization of the pasta 
phases? What property can be used to signal a change 
of pasta phase? The purpose of the present work is to 
construct the instruments needed to properly quantify 
the pasta structures. 

Taking advantage of the microscopic details produced 
by a classical molecular dynamics model, this investiga- 
tion combines the power of cluster detection algorithms 
used in nuclear collisions with indicators borrowed from 
condensed matter physics and topology to detect the 
transitions between pasta phase structures in a quanti- 
tative way. 

After an introduction of the model, we will introduce 
a series of techniques used to classify the pasta struc- 
tures that will help us reach our goal. Starting from 
global measures to understand the cluster composition 
(fragment size distribution, nucleon mobility and persis- 
tence, fragment isotopic composition and radial distri- 
bution function), we will progress into topological tools 
(Minkowski functionals) that will allow us to character- 
ize the shape of the pasta structures as well as to detect 
changes between them. A final discussion of the results 
will then help us reach a series of conclusions and to draw 
an outlook of the future tasks. 



II. NUCLEON DYNAMICS 

To study the structure of stellar crusts is necessary 
understand the behavior of nuclcons at the proper den- 
sities, temperatures and proton-to-neutron ratios; such 
knowledge comes from the study of heavy ion fragmenta- 
tions. The initial statistical studies of nuclear collisions 
of the 1980's [29|,[3(J, rapidly gave way to dynamical the- 
ories based on classical, semiclassical and quantum ap- 
proaches. 

The semiclassical models use the Boltzmann-Uehling- 
Uhlenbeck equations [3l| to track the time evolution in 
phase space of the probability of finding a particle mov- 
ing in a mean field. On the other hand, the quantum 
molecular dynamics models (QAID) solve the equations 
of motion of nucleon wavepackets moving within mean 
fields. Unfortunately, these theories either do not lead 
to cluster formation or yield a poor description of cluster 
properties and both must resort to the use of all sorts 
of extraneous techniques such as adding fragments by 
hand, coupling to "afterburners" to produce secondary 
decays [32|, and introducing hidden adjustable parame- 
ters (e.g. width of wavepackets, number of test particles, 
modifications of mean fields, effective masses and cross 
sections, etc.) to satisfy the operator's taste. 

These problems are either non-existent or much re- 
duced in classical models. Classical dynamical models 
generally solve Newton's equations of motion to track 
individual nucleons moving under two-body potentials; 
coupled to cluster recognition algorithms these calcula- 
tions yield microscopic views of nuclear reactions as well 
as of nuclear structures. The only apparent disadvantage 



of the classical models would be the lack of quantum 
effects, such as the Pauli blocking; fortunately, in stel- 
lar environments the very small nucleon energies lead to 
frozen-like structures where the blocking of momentum- 
transferring collisions ceases to be relevant. [See Ref. [33| 
for a calculation of the mean thermal wavelength of a nu- 
cleus in stellar conditions to justify the use of a classical 
approach.] 

Let this rather long preamble serve to justify extending 
the use of a classical model designed for nuclear reactions 
to the study neutron star crusts. In this work we use a 
classical model to obtain a detailed microscopic picture 
of the pasta structures and be able to detect transitions 
between phases. 



A. Classical Molecular Dynamics 

We use a molecular dynamics code combined with al- 
gorithms for cluster reco gnit ion. Our classical molecular 
dynamics model, CMD [341 . is based on the pioneering 
work of Pandharipande [35[ and has been very fruitful in 
nuclear studies of, among other phenomena, neck frag- 
mentati on 13611 . phase transitions 1371. 13811 , critical phe- 
nomena [33 . |40( , the caloric curve [4ll . |42| , and isoscaling 
(iH |44| all without any adjustable parameters. Readers 
are directed to these references for further details on the 
model; here only a brief synopsis will be presented along 
with its extension to infinite systems. 

In a nutshell, CMD treats nucleons as classical parti- 
cles interacting through a two-body potential and solves 
the coupled equations of motion of the many-body sys- 
tem to obtain the time evolution of all particles. Since 
the (r, p) information is known for all particles at all 
times, it is possible to know the structure of the nuclear 
medium from a microscopic point of view. 

CMD uses the phenomcnological potentials developed 
by Pandharipande [35| : 

V np (r) = V r [exp(-fx r r)/r - exp(-fx r rc)/r c ] 

- V a [exp(~p a r)/r - exp(-p a r a )/r a ] 
V N N{r) = Vo[exp(—fjt r)/r — exp(-fior c )/rc] , 

where V np is the potential between a neutron and a pro- 
ton and it is attractive at large distances and repulsive at 
small ones, and Vjvjv is the interaction between identical 
nucleons and it is purely repulsive. Notice that no bound 
state of identical nucleons can exist, also notice that, at 
a difference from potentials used by other models [la ], 
these potentials have a hard core. 

The cutoff radius is r c = 5.4 fm after which the poten- 
tials are set to zero. Two sets of values for the Yukawa 
parameters /Lt r , fj. a and po were fixed by Pandariphandc 
to correspond to infinite-nuclear matter systems with an 
equilibrium density of po = 0.16 /m -3 , a binding energy 
E(po) = 16 MeV/nucleon and compressibility of about 
250 MeV ("Medium") or 535 MeV ("Stiff") 
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the past, a combination of Monte Carlo and molecular 
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FIG. 1: (Color online) Smorgasboard of pasta shapes cor- 
responding to the densities shown and to x = 0.5 and 
T = 0.1 MeV. 

dynamics techniques was applied within a statistical for- 
malism to obtain neutron start crust properties (TEj . 

B. Simulating the neutron star crust 

To study the neutron start crust we use CMD to sim- 
ulate an infinite medium. Systems with 2000 or 3000 
nucleons were constructed and replicated with periodic 
boundary conditions in 26 surrounding cells. In particu- 
lar, the proton ratios used were of a: = Z/A = 0.5 (f 000 
neutrons and 1000 protons) or 0.3 (2000 neutrons and 
fOOO protons). The cubical box size used was adjusted 
as to achieve densities between p = 0.01 /m~ 3 (about 
Po/15) and p . 

As the crust is expected to be embedded in a degen- 
erate electron gas produced by weak decays during the 
supernova explosion, it is necessary to take into account 
its Coulomb interaction. Although the nucleon-electron 
system is overall neutral and /3-equilibrated, the infinite 
Coulomb range requires the use of some approximation; 
two common approaches are the Thomas-Fermi screened 
Coulomb potential (used in QMD in [l7j) or the Ewald 
summation procedure [45[. Although CMD is able to op- 
erate under either approximation, in this work the former 
is adopted (see [46[ for a comparison of methods under 
CMD). 

Approximating the electron g uniform ideal 

Fermi gas at the same number density as the protons, 
its effect can be included in the nucleon's equations of 
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FIG. 2: Temperature evolution of the fragment size distri- 
bution obtained from 200 configurations with x = 0.3 and 
p~po/10. 



motion by means of the screened Coulomb potential ob- 
tained from the Poisson equation: 

V# Cr) (r)^cxp(-r/A) 

where the relativistic Thomas-Fermi screening length is: 

_ j_ 

A = (7r 2 /2e) \kp (k'p + ml)] 4 , m e is the electron mass, 

1/3 

kp = {Zir 2 p e ) is the electron Fermi momentum, and 
p e is the electron gas number density equal to that of the 

1 /3 

protons. The size of the simulation cell, L = (A/p) 1 , 
should be significantly larger than A; in our case we sat- 
isfy such requirement using the prescription of [l5[ and 
setting A = 10 fm. 

The trajectories of individual nucleons, now governed 
by the Pandharipandc and the screened Coulomb po- 
tentials, are then tracked using a Verlet algorithm with 
energy conservation of 0(0.01%). The system is force- 
heated or cooled using isothermal molecular dynamics 
with the Andersen thermostat procedure [47[ which grad- 
ually cools in small temperature steps while reaching 
thermal equilibrium at every step. We focus in the range 
of T = 0.1 to 1.0 MeV] although this last temperature is 
large for stellar crusts, in terms of the nucleon dynamics 
it practically corresponds to a frozen state. 



III. CHARACTERIZING THE CRUST 

At a difference from most QMD simulations, which 
tend to track individual evolutions, here we obtain reli- 
able statistics by sampling 200 times each configuration 
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Relative abundance of large clusters as a function of 
for x = 0.3 and 0.5 at T = 0.1 MeV. 
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FIG. 4: Root mean square displacement of the nucleons from 
their original position as a function of the simulation time 



steps for systems with p 
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0.015 /m" d and x = 0.3 at the 



with specific x, p, and T conditions. Figure [T] shows 
an smorgasboard of Italian delicacies produced by CMD 
with x = 0.5, T = 0.1 MeV and twenty different densi- 
ties; please notice that for clarity the figures do not show 
single nucleons, i.e. the gaseous phase. In spite of their 
beauty, one cannot use those figures to properly charac- 
terize the pasta shapes, for that one must resort to other, 
less visually attractive, techniques. 

On each of the configurations achieved, the nucleon 
(r, p) information is recorded and used later to identify 
clusters and to characterize the structure by means of the 
liquid structure function and the Minkowski functionals. 



A. Cluster composition 

The nucleon positions are used to identify clusters by 
means of the "Minimum Spanning Tree" (MST) algo- 
rithm refined for nucleon dynamics in j48l l49j . In sum- 
mary, MST looks for correlations in configuration space: 
a particle i belongs to a cluster C if there is another par- 
ticle j that belongs to C and | r j — r j\ < r c i , where r c i is 
a clusterization radius which is set to 3.0 fm. In spite 
of using only r-space correlations, MST yields accurate 
results in the case of stellar crusts due to the low tem- 
peratures and small momentum transfer, and thus here 
it is preferred over other more robust cluster-detection 
algorithms (such as the "Early Cluster Recognition Al- 
gorithm", ECRA (5t| , which take into account relative 
momenta and binding energies). In our case of periodic 
boundary conditions, the MST method was modified to 
recognize fragments that extend into adjacent cells. 

Figure [2] shows examples of the fragment population 
obtained with CMD -MST at x = 0.3, p = 0.015 fm' 3 



and four different temperatures; of particular interest is 
the evolution of the clustering as a function of the tem- 
perature as it shows a trend opposite to that observed in 
nuclear collisions. The figure shows a typical evolution 
for the range 0.1 MeV < T < O.lMeV and, as it can 
be clearly observed, the large fragment multiplicity in- 
creases with T. This is in opposition to what happens in 
heavy ion reactions at high energies where heavier frag- 
ments shrink in size by particle evaporation during the 
final expansion stage of the reaction. In the case of in- 
finite systems, however, the lack of expansion (and lack 
of a reduced pressure) makes evaporation less probable 
and, combined with the possibility of connecting frag- 
ments to neighboring cells, it favors the growth of cluster 
sizes as soon as the nucleons reach enough mobility with 
increasing T. 

The growth of large fragment multiplicity can also be 
seen as a function of the density. Figure [3] shows a typi- 
cal behavior of the relative multiplicity of large clusters, 
A > 100 and A > 10, obtained at different densities for 
both x = 0.3 and 0.5 at T = 0.1 MeV, As the density 
increases, the number of clusters of A > 100 increases 
practically linearly with p up until a single large frag- 
ment is formed; smaller clusters of A > 10 are abundant 
at low densities but decrease for larger densities. The 
density at which the number of large clusters condense 
into a single one can be thought of as a "percolation" 
density, this value, of course, depends on the simulation 
parameters such as number of particles, cell size, temper- 
ature, etc.; for the cases shown, the percolation densities 
are p « 0.03 fm~ 3 for x = 0.3 and p sa 0.024 fm~ 3 for 
x = 0.5 at T = 0.1 MeV. 

The dynamics of the nucleons within systems in equi- 
librium can be gauged through their average displace- 
ment as a function of "time" , i.e. through the time steps 
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T=0.1 MeV 




FIG. 5: Time evolution of the persistency for the listed tem- 
peratures as a function of the simulation time steps for x — 0.3 
and p = 0.015 fm~ 3 . 



of the simulation. Figure [4] shows the root mean square 
displacement of the nucleons from their original position 
in 200 time steps; as a metric one must remember that 
the range of the potential is 5.4 fm and -as we will see in 
section IIII Bl - the interparticlc distance at these densi- 
ties is of the order of 1.7 fm. The increment in mobility 
as a function of the temperature is obvious. 

Likewise, the microscopic stability of the clusters can 
be quantified through the "persistency" [5l|, [52| which 
measures the tendency of members of a given cluster to 
remain in the same cluster. Figure [S] shows the time evo- 
lution of the persistency for systems with p = 0.015 fm~ 3 
and x = 0.3 at the listed temperatures; notice that a per- 
sistency of ~ 1 indicates that most of the particles remain 
in the same cluster, while smaller values indicate a larger 
exchange rate. The anti-correlation between this and the 
previous figure is clear, more mobility implies less persis- 
tency, and viceversa. 

Another interesting descriptor is the isospin content of 
the clusters produced. By keeping track of the number 
of protons and neutrons on each fragment it is possible 
to determine the x value for each cluster. An example of 
this is shown in figure |6] where the x content of the frag- 
ments is plotted as a function of the mass of the clusters 
obtained at a density of p = 0.015 fm~ 3 and with x = 0.5 
(top two plots) and x = 0.3 (bottom four). Several ef- 
fects are noticeable: small clusters (A < 10) tend to have 
less protons than the average resulting in smaller x val- 
ues; for the case of a; = 0.3 there is a prominent binding 
of one proton to two neutrons which results in clusters 
of all sizes with values of a; w 1/3; the previous effect is 
not present in the case of x = 0.5 in which all the cluster 
maintain their x values around 1/2. 

We close this subsection noticing that, in spite of be- 
ing a good indicator of the percolating density, the clus- 
ter multiplicity is a poor descriptor of the pasta shapes. 
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FIG. 6: Example of the x value (a; = p/A) of the clusters ver- 
sus their size, corresponding to systems with p = 0.015 fm~ 3 
and x = 0.5 (top two plots) and x = 0.3 (bottom four plots) 
at the temperatures listed. 



Stepping up in complexity, we now turn to the radial 
correlation function to probe the pasta al dente. 



B. Pair correlation function 

Another global characterization of the structure of nu- 
clear matter is obtained from the pair correlation func- 
tion, g(r), which is the ratio of the average local density 
to the global density, g(r) = p(r)/po; it gives information 
about the spatial ordering of the nuclear medium. 

For computing purposes, the pair correlation function 
g(r) is taken as the conditional probability density of 
finding a particle at + r given that there is one particle 
at r,. Formally, 



9(r) 



V 



4irr 2 N 2 




where r,j is |r — For our case, this was calculated by 
constructing histograms of the distances between parti- 
cles for several configurations obtained with the same x, 
p and T and then averaging them; to include all particles 
and their images the range was extended to < 1.5L. 
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FIG. 7: Examples of the radial correlation function for x = 
0.5, T = 0.1 MeV and various densities. Also shown is the 
case when nearest neighbors are as probable as second nearest 
neighbors indicating a "lasagna" type structure; r is in fm. 



Figure [7] shows examples of g(r) obtained for cases 
with 2,000 nucleons, x = 0.5, T = 0.1 MeV and den- 
sities 0.01 < p < 0.039 fmT z . The inset shows the case 
p = 0.048 /to -3 when the nearest neighbors are just as 
probable as second nearest neighbors signalling the onset 
of a "lasagna" type structure, cf. Fig. [TJ 

It is worth noticing that, in the case shown, the lo- 
cation of the nearest neighbors remains constant at r ~ 
1.7 fm at all densities. This is due to the fact that -since 
at subcritical densities the medium is metastable or un- 
stable, it breaks into a gaseous and a condensed phasc- 
the position of the peak of g(r) is an average between the 
location of neighbors in the gas-liquid mixture; the con- 
densed matter at normal density has nearest neighbors 
at r ps 1.4 fm [Hj. 

Once again, in spite of the rich information derived 
from g(r), it is still insufficient to tag the phases unequiv- 
ocally; for this, other more complex constructs must be 
borrowed from cosmology and, ultimately, from topology. 



C. Topological constructs 

The most obvious properties of closed surfaces that can 
be used to characterize their shapes are the volume V, 
surface area A, and the curvature. The latter is less triv- 
ial than the other two as it does not hold a unique value 
for a given shape but varies from point to point; however, 
mean curvatures of a closed body can be obtained by an 
averaging procedures such as the "integral mean curva- 
ture", defined as H = J df {Ri + R 2 )/2R 1 R 2 where R x 
and i?2 are the principal radii of curvature of the surface 
and df is a differential of area. 

In general, V, A, H plus an interesting construct 
known as the Euler characteristic are collectively known 
as the "Minkowsky functionals" ; according to integral ge- 



FIG. 8: (Color online) Sample transformation of a nuclear 
structure to a corresponding polyhedron. The structure cor- 
responds to a case with x = 0.5, p = 0.33 fm~ 3 and 
T = 0.1 MeV. 
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FIG. 9: (Color online) Spatial configurations formed under 
T = 0.4 and p = 0.045 /m~ 3 (left) and T = 1.0 and p = 
0.072 fm~ 3 (right), both for x = 0.3. 



ometry the morphological properties of ZD objects can 
be completely described in terms of them. 

The XVIII century work of Euler-L'Huilier showed 
that, independent of the shape of any polyhedra, when 
deducting the number of edges from the number of ver- 
tices and adding the number of faces it always yields 2 
plus twice the number of cavities, quantity now known 
as the "Euler characteristic" , x- Although this previous 
property is for solids bounded by plane surfaces, it also 
applies in any 3D surfaces with x related to the total 
curvature of the surface through the Gauss-Bonnet the- 
orem. 

In topological terms, two orientable closed surfaces are 
homcomorphic to each other if their Euler characteris- 
tics are the same. Conversely, two homomcorphic closed 
surfaces will always have the same value of x- There- 
fore, since our pasta niblets are all orientable, their Euler 
characteristic completely classifies them up to an homeo- 
morphism; adding the rest of the Minkowsky functionals 
eliminates such redundancy and guarantees a complete 
classification of the pasta shapes. 

An immediate problem is the fact that the nuclear clus- 
ters are not polyhedra and do not even form closed sur- 
faces. This obstacle, however, can be circumvented by 
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FIG. 10: Multiplicites corresponding to figure|9j although not 
identical they are very similar. 



replacing the nuclear structure with a scaffold-like arma- 
ture composed of cells enclosing each one nucleon. In our 
case this is done through the algorithm of Michielsen and 
De Raedt [HH which has already been used in the study 
of stellar crusts albeit in a different methodology [2l| . 

Synoptically, the simulation volume is subdivided into 
a mall of cubic cells. Those cells which contain the coor- 
dinates of a nucleon are kept while the rest are deleted. 
The sizes of the cells are made smaller than the nearest 
neighbor distance found in g(r) to enforce a one-particlc- 
per-cell occupation, but not too small as to avoid creating 
spurious cavities between neighboring nucleons. It is on 
this imaginary platform that the Minkowsky functionals 
of the nuclear structure are computed. 

In general, x equals the number of regions of connected 
grid cells minus the number of completely enclosed re- 
gions of empty grid cells. Two grid cells are connected 
if they arc immediate neighbors, next-nearest neighbors, 
or are connected by a chain of occupied grid cells. Char- 
acterizing the connected structure by its number of oc- 
cupied cubes, n c , edges, n e , faces, n/, and vertices, n v , 
including possible contributions from the interior of the 
structure, the Minkowski functionals can be calculated 
through [53| 



IB = 3n c 



V = n c , A 
2n,f + n c , x z 



-6n c + 2n/, 
-n c + Hf — n e 



where V stands for the volume, A for the area, B for the 
mean breadth B and \ f° r the Eulcr number; the mean 
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FIG. 11: Euler and curvature of the structures shown in fig- 
ure [9] the difference between the two cases is easy to spot. 



breadth measures of the average "size" of a body and it 
is related to the integral mean curvature H mentioned 
before. Figure |8] shows a typical nuclear structure along 
with the grid constructed around it; the values of the 
Minkowski functionals obtained from such grid are Cur- 
vature = 215 and Eulcr = —17, as we will see next such 
shape can be classified as a "jungle gym" . 



D. Topological classification of the pasta 

To illustrate the use of topology to classify the pasta 
shapes let us use two seemingly similar structures ob- 
tained with x = 0.3 but at different densities and tem- 
peratures, namely p = 0.045 /to" 3 and T = 0.4 MeV 
and p = 0.072 /m" 3 T = 1.0 MeV. The spatial config- 
urations of these two cases, practically identical to the 
eye, are presented in figure Although there are minor 
differences (see insets), figure [TD] shows that both config- 
urations have very similar mass multiplicities. 

The difference between the configurations, however, 
surfaces when wc calculate the corresponding curvature 
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FIG. 12: (Color online) Typical artificial "pastas" used to test 
the classification powers of the Euler- Curvature combination. 
Respectively they are the "jungle gym" (top left), "lasagnas" 
(top right), "straight spaghettis" (bottom left) and a "curled 
spaghetti" . 



and Eulcr numbers. Figure [Til shows the evolution of the 
curvature and of the Euler number as the simulations ad- 
vances after equilibration. Clearly shown are the inverted 
values of Euler and the curvature for the two cases; while 
the curvature is positive (with less cavities than bellies) 
and large in the low T low p case, it becomes negative 
(with more internal cavities than tummies) and smaller 
in the opposite case. 

To investigate this point further, we created artifi- 
cial structures in the form of gnocchi, spaghetti, lasagna 
and crossed-lasagnas, which we call "jungle gym" , and 
their inverse structures (with voids replacing particles 
and viceversa) and calculated the values of the two topo- 
logical variables; some of the structures used are shown in 
figure [T2] and their location in the Curvature- Euler plane 
in figure Q2] The magnitudes shown are determined by 
the size of the structure as well as the digitization cell 
size. In general one can conclude that lasagnas tend to 
lie near the origin, spaghettis have near zero Euler num- 
bers and positive curvatures, gnocchis have positive cur- 
vatures and Euler numbers, and "jungle gyms" positive 
Euler number and negative curvature; all anti-structures 
reverse the curvature but maintain the Euler character- 
istic. All cases calculated at all x values, densities and 
temperatures were observed to satisfy this classification. 

For instance, the structure in figure [8] with curvature 
215 and Euler number —17 is clearly a "jungle gym"; it 
must be remarked that -to our knowledge- this is the 
first time this type of structure has been reported. Like- 
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wise, the structure in the left panel of figure [^1 with both 
negative curvature and Euler number, can be classified as 
an "anti-jungle gym" , whereas the accompanying struc- 
ture with positive curvature and negative Eulcr number 
would have to be classified as a "jungle gym" . 

Table UJ shows the classification of several of the struc- 
tures obtained in our study at T — 0.1 MeV and at the 
listed densities, and figure ITU shows their location in the 
C — E plane; the x = 0.5 column of the table and the 
circular dots on the figure correspond to the structures 
presented in figure [TJ As the absolute magnitude of the 
curvature and Euler number depends on the overall size 
of the structure, i.e. on the number of particles used, 
the data in table [X] were normalized to have maximum 
absolute values of 100. In the table, the classifications 
are abbreviated as G for gnocchi, J for jungle gym, L for 
lasagna, S for spaghetti, and AG, AJ, AL and AS for the 
reverse structures. 



IV. CONCLUDING REMARKS 

And thus we have reached our objective. The com- 
bination of curvature and Eulcr number and a proper 
recognition of fragments appears to be robust enough as 
to uniquely classify the shapes attained by the nucleons 
at densities, temperatures and isospin content of inter- 
est in the study of neutron star crusts. In obtaining this 
result, the classical molecular dynamics and associated 
tools (cluster recognition algorithms, persistence, etc.) 
proved to be a valuable method, which we now plan to 
exploit. 
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In future investigations we will apply this method to 
an in-depth study of the origin of clustering. As hinted 
in [HI, the role of the long range Coulomb interaction 
has not been fully explored; current exploratory runs are 
indicating that Coulomb merely shifts the space scales a 
bit but is not an "if-and-only-if" requirement for the for- 
mation of clusters. Along the same lines, as the present 
CMD model can function both with a stiff and medium 
compressibility potentials, we also plan to investigate the 
role of the equations of state in the formation of the pasta 
structures. 



o 

Curvature 



FIG. 13: (Color online) Systematic classification of the arti- 
ficial structures in terms of the curvature and Euler number. 
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FIG. 14: Curvature-Euler coordinates of the structures listed 
in table UJ Labels "A", "a", "L" and "t" correspond, respec- 
tively, to the initial and final points in the table. 
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